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Abstract. A classical realization of the two-site Bose-Hubbard Hamiltonian, based 
on light transport in engineered optical waveguide lattices, is theoretically proposed. 
The optical lattice enables a direct visualization of the Bose-Hubbard dynamics in Fock 
space. 
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The Bose-Hubbard Hamiltonian provides a paradigmatic theoretical framework to 
investigate the physics of strongly interacting bosonic systems pp. In particular, the 
two-site Bose-Hubbard model has received a great deal of attention in the last decade 
as a simple model to investigate the dynamics of ultracold bosonic atoms confined in a 
double- well potential El SI El El El IB] , the so-called bosonic junction (for a recent review 
see [9]). A remarkable property of the bosonic junction is the existence of Josephson-like 
oscillations of the atomic populations and of macroscopic self trapping due to atom-atom 
interaction, which were originally predicted within a semiclassical (mean-field) limit of 
the Bose-Hubbard model [21 [3]. However, various quantum features of the Bose-Hubbard 
dynamics, such as quantum fluctuations, fragmentation [TU], the cat-states formation 
in the supercritical attractive case [11], or the many-body coherent control of tunneling 
[L2"] [T3"] . are not accessible within the mean-field limit. This has lead to intensive 
studies on the relation between the full quantum and mean-field dynamics (see, e.g., 
[H I [T5l [T6l [T7] and references therein). A comprehensive description of the entire variety 
of phenomena rooted in the Bose-Hubbard Hamiltonian would benefit from a direct 
access to quantum dynamical evolution in the Fock space. However, in experiments 
with ultracold gases the measurable quantities are usually atomic population imbalance 
and relative phases [3], whereas full information on Fock state occupation evolution is 
generally not accessible. 

In another physical context, light transport in waveguide lattices has provided a 
test bench to visualize in photonic systems the classical analogues of a wide variety of 
coherent single-particle quantum phenomena generally encountered in condensed matter 
or matter wave systems [HI [191 EH], such as Bloch oscillations and Zener tunneling 
[TSJ [5TJ , dynamic localization f2~2~\ |2"3"] , and Anderson localization [2"H 123] to mention a 
few. As compared to quantum systems, photonic lattices enable a direct and complete 
visualization of the dynamics by mapping the temporal evolution of the quantum 
system into spatial propagation of light waves in the engineered lattice [19J. Most 
of such previous quantum-optical analogue studies, however, focused on single-particle 
effects, whereas the possibility to visualize in a classical optical setting the analogues 
of many-particle phenomena, such as those rooted in the Bose-Hubbard model, remains 
unexplored. Recently, photonic lattices with engineered coupling constants have been 
proposed as classical analogs to quantum coherent and displaced Fock states [26] . 

In this Letter it is shown that photonic lattices with suitably engineered coupling 
and propagation constants provide a simple realization in the Fock space of two-site 
Bose-Hubbard Hamiltonian. 

The starting point of the analysis is provided by a standard two-site Bose-Hubbard 
Hamiltonian that describes a system of N interacting bosons occupying two weakly- 
coupled lowest states of a symmetric double-well potential, which reads (see, for instance, 



where di^ (a\ 2 ) are the bosonic particle annihilation (creation) operators for the modes 




(i) 
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in the two wells, J > accounts for the coupling constant between the two modes and 
U is the strength of the on-site interaction (U > for a repulsive interaction). If we 
expand the vector state of the system \ip(t)) on the basis of Fock states with constant 
particle number N, i.e. after setting 

|1K*)> = T , = a\ l at N - l \0) (2) 

the evolution equations for the amplitude probabilities to find / particles in the left 
well and the other (N — I) particles in the right well, as obtained from the Schrodinger 
equation ihd t \ip(t)) = H\ip(t)), read explicitly 

i^T = -(«jCi+i + «i_ic { _i) +V1C1 (3) 
(Z = 0, 1, 2, ...,N), where we have set 

k, = Jv/(/ + l)(iV-/) , I/, = | + (AT - J)* - N ] . (4) 

The normalization condition YliLn l c «(^)| 2 = 1 * s assumed. 

The evolution equations for the Fock state amplitudes q can be viewed as formally 
analogous to the coupled-mode equations describing the transport of light waves in 
a tight-binding array composed by (N + 1) waveguides with engineered propagation 
constant shift Vi and coupling rates K[ between adjacent waveguides, in which the 
temporal evolution of the Fock-state amplitudes of Bose-Hubbard Hamiltonian is 
mapped into the spatial evolution of the modal field amplitudes of light waves in the 
various waveguides along the axial direction z (see, for instance, [JSH20])- The fractional 
light power distribution |q| 2 in the various waveguides of the array at the propagation 
distance z thus reproduces the occupation probabilities of the bosons between the two 
sites of the double- well. 

In the optical context, the tight-binding model (3) can be derived using a 
variational procedure starting from the paraxial and scalar wave equation for the electric 
field amplitude <f>(x, z) describing the propagation of monochromatic light waves at 
wavelength A in an array of (N + 1) waveguides with refractive index profile n(x) and 
substrate index n s 

A 2 

i\d z (f) = -—d 2 x (j) + [n s - n(x))(f>, (5) 

Z,TL S 

where A = A/(27r) is the reduced wavelength of photons (for details see, for instance, 
[29J). In particular, to independently engineer the coupling constants K\ and propagation 
constant shifts Vi of the waveguides, one can assume a chain of waveguides with 
equal normalized refractive index profile g(x), but with different index contrasts An; 
(/ = 0, 1, 2, N) and spacing di = xi — xi-i (I = 1, 2, N), i.e. 

N 

n(x) - n s = - 22 ^nig(x - xi). (6) 

1=0 
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Figure 1. (color online) (a) Distribution of waveguide distance (di is the distance 
between waveguides (I — 1) and I, I — 1,2,. ..,9) in a waveguide array composed by 
N + 1 = 10 waveguides that realizes the distribution of the coupling constants defined 
by Eq.(4) for J = 0.0781 mm -1 . The panels (b-d) show the refractive index profile 
n(x) — n s of the waveguide arrays that realize the Bose- Hubbard Hamiltonian for 
increasing values of the interaction strength U: in (b) U — 0, in (c) U — 0.0174 mm -1 , 
whereas in (d) U = 0.1043 mm -1 . 



As discussed in the example below, to realize the Bose-Hubbard lattice (3) the values of 
d\ and An; slightly vary around some mean values d r and An that define a uniform array. 
The waveguide separation di mainly determines the value of the coupling rate ki-±, with 
a characteristic exponential dependence of from di, whereas the index change An; 
mainly defines the propagation constant mismatch V\. It is worth noticing that, in the 
absence of particle interaction, i.e. for U — 0, the lattice model (3) associated to the 
Bose-Hubbard Hamiltonian was previously introduced in the photonic context to realize 
exact spatial beam self-imaging in finite waveguide arrays [27] and shown to belong to a 
rather general class of exactly-solvable self-imaging tight-binding lattices with equally- 
spaced energy levels [28] . A nonvanishing value of the interaction strength U breaks the 
periodic self-imaging property of the array and enables to visualize with light waves the 
rich dynamical features of the two-site Bose-Hubbard Hamiltonian directly in the Fock 
space (see, for instance, [21EHHE1E1EIIH1ISI HH HSJ EH E] and references therein) . Here 
two important dynamical features will be considered for the sake of example, namely 
(i) the transition from Josephson-like oscillations to self-trapping as the interaction U is 
increased [2], [31 E], and (ii) the transition from single-atom to correlated pair tunneling 
of two bosons in a double well potential [H [30] . 
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Figure 2. (color online) Light intensity evolution (snapshot of 



the 



three waveguide arrays of Figs.l(b-d) corresponding to excitation of the left boundary 
waveguide, as obtained by numerical simulations of the paraxial wave equation (5). The 
images in (a), (b) and (c) correspond to the Bosc-Hubbard Hamiltonian for increasing 
values of the interaction strength U (from left to right). 



To visualize the first effect, we specifically design three arrays comprising N + 1 = 10 
waveguides, which mimic the evolution of N = 9 atoms in a double well, for three 
increasing values of the interaction U . The arrays were designed to operate at A = 633 
nm in a transparent glass with a bulk refractive index n s = 1.45 and with a normalized 
waveguide channel profile g(x) = {erf[(x + w)/D x ] — erf [(re — w) / D x ]} /[2erf(w / D x )] 
with channel width 2w = 4 /im and diffusion length D x = 0.3 //m -1 . The distances 
between waveguides di were designed to realize the coupling rates Ki given Eq.(4) with 
J = 0.0781 mm -1 . To determine the distribution of distances di, a reference value 
An = 2 x 10~ 3 of refractive index contrast was assumed, and correspondingly the 
coupling rate k between two adjacent waveguides versus distance d was computed, 
yielding to a good approximation the exponential dependence tz(d) = Kq exp[— j(d — d r )] 
for distances not too far from the reference distance d r = 8 um, where k,q ~ 0.3907 mm -1 
and 7 ~ 0.6 urn -1 . The resulting distance distribution, depicted in Fig.l(a), shows 
that the waveguide distances slightly vary around the reference value d r indeed. By 
modulating the index contrasts Arii of the waveguides around the reference value An, 
three different arrays were then designed to realize the three different interaction regimes 
U = 0, U = 0.0174 mm" 1 , and U = 0.1043 mm" 1 , as shown in Figs.l(b-d) [HI]. Such 
arrays could be fabricated in fused silica by the recently developed femtosecond laser 
writing technique [20J, in which the different refractive index contrasts are obtained 
by varying the speed of the writing laser beam. Figures 2 show the evolution of light 
intensity \(f>(x, z)\ 2 along the three arrays, as obtained by numerical simulations of Eq.(5) 
using a standard pseudospectral split-step method, when the left boundary waveguide 
is excited at the input plane. Such an excitation corresponds to the initial conditions 
q(0) = 5i t o, i.e. to the entire iV bosons in the right well at the initial time. Note that 
in case U = [Fig. 2 (a)] periodic self- imaging of the light pattern is observed because of 
the equispacing of the energy levels of the Bose-Hubbard Hamiltonian. Such periodic 
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Figure 3. (color online) Behavior of the numerically-computed population imbalance 
P{z) for the three waveguide arrays of Fig. 2 (solid curves), and corresponding behavior 
predicted by the tight-binding Bose-Hubbard Hamiltonian (3). From (a) to (c), the 
interaction strength takes the values U = 0, U = 0.0174 mm -1 , and U = 0.1043 mm -1 . 

oscillations of the light intensity pattern, previously predicted in Ref.[27j and referred 
to as harmonic oscillations, are thus the optical analogue of the bosonic Josephson 
oscillations in the non- interacting regime [2J. A quantitative measure of the Josephson 
oscillations, generally used in the atom optics context, is provided by the normalized 
population imbalance P(z), which is defined as the normalized difference between the 
average numbers of atoms in the two wells, i.e. P(z) = ^2i =0 [(N — 2l)/N]\ci\ 2 . Note 
that in an optical experiment P(z) can be retrieved from a measurement of the beam 
center of mass position (x(z)) = J dxx\4>(x, z)\ 2 / J dx\<f>(x, z)\ 2 using the simple relation 
P(z) — 1 — 2(x(z))/(Nd r ). For the U = case, the evolution of P(z) along the 
array is depicted in Fig. 3 (a) and compared to the exact curve obtained from the tight- 
binding lattice model (3). Note that P(z) oscillates between -1 and 1 with spatial 
period zr = n / J ~ 4 cm, which is precisely the period of Josephson oscillations of the 
bosonic junction in the absence of interaction. As the interaction strength is increased, 
the self-imaging property of the array is clearly broken [see Figs. 2(b) and (c)], with 
a clear transition to damped Josephson oscillations [Fig.3(b)] and to the self-trapping 
regime [Fig.3(c)]. It should be noted that the mean-field limit of the two-site Bose- 
Hubbard model in the large N limit, which is described by two nonlinear coupled mode 
equations j2J |9j, could be realized in an optical setting by a simple nonlinear optical 
directional coupler [32J, and self-trapping phenomena in such nonlinear couplers have 
been previously investigated. However, our waveguide lattice system provides an exact 
realization of the Bose-Hubbard model even for a low number of particles (as in the 
example previously discussed) where the mean field approximation fails. 

As a second example, we discuss the phenomenon of correlated pair tunneling of 
two bosons in a double well potential induced by interaction, which was investigated 
theoretically and experimentally in previous works (see, for instance, [HI [30]). Even if the 
two-mode Bose-Hubbard model is capable of describing the dynamics of such a system 
solely for weak interactions and far from the fragmentation regime [30], it well explains 
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the transition from Rabi oscillations of single particles to correlated pair tunneling in 
the presence of a weak interaction. Such a kind of correlated dynamics of a pair of 
interacting atoms was experimentally observed in Ref. [8] by recording both the atom 
position and phase coherence over time, and was explained on the basis of the simplified 
Bose- Hubbard Hamiltonian (2). In the N = 2 particle case, the coupled-mode equations 
(3) reduce to the following ones 
.dc 



' dt 
.dc\ 

'~dt 
.dc 2 
'~dt 



- V2Jd + Uc 

- V2J(c Q + c 2 ) 

- V2Jci + Uc 2 . 



(7) 



The tunneling dynamics can be analyzed by the introduction of the percentage of bosons 
in the right well, puit) = | c ) | 2 + ( 1/ 2) | ci (t) | 2 , and the pair (or same-site) boson 
probability P2(t) = \co{t)\ 2 + | C2 ) | 2 [30], i.e. the probability to find the two bosons in 
the same well (either the left or the right well). The evolution of the particle occupation 
probabilities |c;(t)| 2 can be readily obtained by solving Eq.(7) with the initial condition 
co(0) = 1 and ci(0) = 02(0) = 0. One obtains 

2J 2 



\ci(t)\ 2 

\02(t)\ 2 



M 2 
1 



sin 2 (Mt) 



U 
M 



1 + cos 2 (Mt) - 2 cos ( — ) cos(Mt) 

sin 2 (Mt) 



sin 



Ut 

Y 



s'm(Mt) + 



Ut\ 

2" J 

U 2 



AM 2 



(8) 



\c (t)\ 2 = 1 - \ Cl (t)\ 2 - \c 2 (t)\' 



where we have set M 
reads 

P R {t) = 1 



a/4J 2 + U 2 /A. Correspondingly, the behavior of pn(t) and p^it) 
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MMt) + -^— sin 2 (Mt) 



AM 2 



P2(t) 



2J 2 

M 2 



sin 2 (Mt) 



(9) 



(10) 



A typical behavior of both pn(t) and p^if) for increasing values of the ratio U/J is 
shown in Figs. 4(a) and (b), respectively. Note that, for the non-interacting case (U = 0, 
curve 1 in Fig. 4) the atoms simply Rabi oscillate back and forth between both wells, 
and they tunnel independently. As a correlation is introduced, the tunneling becomes 
a two-mode process and the tunneling period increases [curves 2 and 3 in Fig.4(a)]. 
More interestingly, from the behavior of p2{t) one can see that as the interaction is 
increased both atoms remain essentially in the same well in the course of tunneling, 
i.e they tunnel as pairs and the amplitude Ci(t) gets negligible [30]. In our optical 
setting, such a dynamical behavior simply describes light tunneling among three coupled 
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Figure 4. (color online) Behavior of (a) the percentage of bosons in the right well 
Pn{t), and (b) of the pair (same-site) boson probability P2(t) in the Bose-Hubbard 
model with two atoms for for increasing values of the intercation strength U, normalized 
to the hopping rate J. Curve 1: U/J = 0; curve 2: U/J — 4; curve 3: U/J — 8. 

waveguides, with the same coupling rate but with the propagation constant of the 
central waveguide detuned from that of the outer waveguides [see Eq.(7)], the detuning 
being proportional to the interaction strength U in the quantum mechanical analogue. 
Such an optical structure was recently proposed and realized for the observation of an 
optical analogue of two-photon Rabi oscillations in Ref . [33] . Indeed, the large interaction 
regime \U/J\ 3> 1 of the Bose-Hubbard model, corresponding to \ci(t)\ <C 1 and to pair 
tunneling (also referred to as second-order tunneling [5]), leads to Rabi oscillations of 
light power between the outer waveguides of the triplet system, with small excitation of 
the central waveguide. 

In conclusion, a photonic realization of the two-site Bose-Hubbard Hamiltonian 
in engineered waveguide lattices has been proposed. Such an optical setting enables 
to visualize in the Fock space the main dynamical aspects of interacting bosons. 
In particular, waveguide lattices have been designed to visualize the transition from 
Josephson-like oscillations to self-trapping, as well as the transition from single-atom 
to correlated pair tunneling in a simple two-boson system. It is envisaged that the 
idea proposed in this work to use engineered waveguide lattices to simulate in a purely 
classical setting the quantum physics of interacting particles should motivate further 
theoretical and experimental studies. For example, longitudinal modulation of the 
refractive index in the lattice or the introduction of balanced loss and gain in the 
waveguides might be used to mimic the recently proposed kicked Bose-Hubbard [T2] 
and non-Hermitian PT-symmetric Bose-Hubbard models [MISS]. 

This work was supported by the Italian MIUR (Grant No. PRIN-20082YCAAK). 
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